Subroutine SIMULATION(EXP_L_prime,B, RET_PRIME, jExperimentNumber )

Use PARAMETERS
USE Functions_libs
USE Interpolate_func
USE random

Implicit NONE

Integer :: jzt_int, jo_lagged, juntreated, joption_pay, jr_rep, jtype2_new, jdp_new, jdu_new, js_new,  I_pay,  jhc_ret, jmc, jtrpay, I_treat, jos, jtype4, jtype4_p, jno, jBenchmark_exp, jh_alt, jr_alt, jr_p_alt, jexp_new, jMEcat_p, jMEcat, jExperimentNumber , jspouse,counter, jzhc_p,jzhc, jmar,jmar_p, jtype3, jtype3_p, jzt, jzt_p,jsample,jas,jas_p,jsav,je,jh,jh_p, jr,jr_p,ji, jage,j, jht, jyears, jfix, jtype1,jtype2,jtype2_p,js,jdp,jdu,jx, jo,jo_p, jins, jhrs , jx_age, jo_act , jsav_act, kk  , jsav_upper , TR_indic, jinv_upper  , jsav_upper_ret  , jhk, jinc


Real (Kind=dd), dimension (n_as,n_type1,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime  !value functions    
Real (Kind=dd), dimension (n_o,n_type4,n_zt,n_as,n_type3,n_type1,n_hc_grid,n_ret-1)::  B !value function while working. (1=not have option to treat and not pay, 2=have option)

Real (Kind=dd), dimension (n_h,n_r,n_o,n_mar,n_zhc)::  EXP_L_prime_INT


Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET !value functions when retired
Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_ret : n_age,n_mar):: RET_PRIME !value functions when retired
Real (Kind=dd), dimension (n_MEcat,n_h,n_r,n_type2,n_mar):: RET_PRIME_INT


Real  (Kind=dd) :: Resources, prob_s, prob_d, income_bt, income, base_1, base_2 , Total_tax, actual_earn, tax_SS, tax_med, V_death, V_c_all_trpay, V_cf_trpay
Real (Kind=dd) :: last_value, maxB, valret   , check
Real (Kind=dd) :: UTI,x, EXPECTED_RET_PRIME , UTI_print, V_death_print

Real (Kind=dd), dimension (1:443)::  ran_hypo ! these are the saved random draws from the simulation without shocks.
Real (Kind=dd), dimension (1:213)::  ran_hypo_old

Real (Kind=dd):: CON, TR, ASP , MAX_SAV  , as_opt, TR_opt ,   INT_sav_upper_ret   , INT_sav_upper, EXPECTED_B_PRIME
Real (Kind=dd), dimension (n_o) :: EXPECTED_B


Real (Kind=dd) :: jx_int_p, prob_sim
Real (Kind=dd), dimension (n_zhc) :: HK_INTERPOL_sim

Real (Kind=dd):: Z_INT_ALT , Z_INT_OFFER_ALT , temp,  jcons_act, jas_p_act, jtr_act  ,jassets  , jx_int, Z_INT , Z_INT_OFFER , B_act , RET_act  , jx_age_interpol  , EXPECTED_RET_PRIME_act, EXPECTED_B_PRIME_act
Integer :: bb,d,o,f   , gr, gg
Real (Kind=dd) :: rmaxipzt  , rmaxipzhc
Real (Kind=dd), dimension (2,n_MEcat,n_type2,n_o) ::  B_int

Real (Kind=dd), dimension(3)::  CON_sim,  ASP_sim, B_sim,   TR_sim
Real (Kind=dd), dimension(2):: income_pay, Total_tax_pay, base_2_pay, base_1_pay
                    
Integer :: indicator_ret, idum, ind_death, intepsi
Real (Kind=dd), dimension(n_mar) :: rmaxiTPzp
Real (Kind=dd) ::  rmaxiIDzp
Real (Kind=dd) :: rx, prob, ry, rb, rd, rf, ri , ro, rinv, i
Real (Kind=dd) :: epsi,weight_epsi, zp_step


1008 Format(I5,1x,I2, 1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,F6.2 , 1x, F6.2 , 1x, F9.2 , 1x,F9.2, 1x,F9.2, 1x,F10.2, 1x,F10.2, 1x, F9.2, 1x, F8.2 , 1x, F10.2, 1x, F4.1, 1x,F9.2, 1x, F10.2,1x,I1,1x, F9.2,1x,F10.2,1x,F10.2,1x,F10.2,1x,F10.2, 1x, I1, 1x, I1, 1x, I1, 1x, I1, 1x, I1 , 1x, F10.6,1x, F10.6,1x, F10.6)
 1005 Format(F16.11) ! this is for random shocks 

!Grid for savings is finer in the simulation     
sav_sim(:)=0.0d0
open(258, file='External_Parameters/savings_grid_sim.txt')
do jsav=1, n_sav_sim
   read (258,*) sav_sim(jsav)
end do 
close (258)
sav_sim(:)=sav_sim(:)/5200.0d0     
  

! sav_grid: useful for interpolation to find the maximum savings
sav_grid_sim(:)=0.0d0
 do jsav=1, n_sav_sim
 sav_grid_sim(jsav) = dble(jsav)
 end do
     
! NUMBER OF PEOPLE OF EACH TYPE SIMULATED
! I simulate histories of n_sample agents. In_type1(jtype1)  tells how many agents of each type are simulated. (There are 18 types in total). 
do jtype1=1, n_type1
n_samp(jtype1)= dble(n_sample)*In_type1(jtype1)  
end do

!***********************************************************************
! Simulating the economy
!***********************************************************************


open(14, file='Simulated_Data/Random_numbers0.txt')
open(15, file='Simulated_Data/Random_numbers_old0.txt')


 if (jExperimentNumber == 1) then
open(12, file='Simulated_Data/Simulated_data1.txt')
else if (jExperimentNumber == 2) then
open(12, file='Simulated_Data/Simulated_data2.txt')
else if (jExperimentNumber == 3) then
open(12, file='Simulated_Data/Simulated_data3.txt')
else if (jExperimentNumber == 4) then
open(12, file='Simulated_Data/Simulated_data4.txt')
else if (jExperimentNumber == 5) then
open(12, file='Simulated_Data/Simulated_data5.txt')
else if (jExperimentNumber == 6) then
open(12, file='Simulated_Data/Simulated_data6.txt')
else if (jExperimentNumber == 7) then
open(12, file='Simulated_Data/Simulated_data7.txt')
else if (jExperimentNumber == 8) then
open(12, file='Simulated_Data/Simulated_data8.txt')
else if (jExperimentNumber == 9) then
open(12, file='Simulated_Data/Simulated_data9.txt')
else if (jExperimentNumber == 10) then
open(12, file='Simulated_Data/Simulated_data10.txt')
else if (jExperimentNumber == 11) then
open(12, file='Simulated_Data/Simulated_data11.txt')
else if (jExperimentNumber == 12) then
open(12, file='Simulated_Data/Simulated_data12.txt')
else if (jExperimentNumber == 13) then
open(12, file='Simulated_Data/Simulated_data13.txt')
else if (jExperimentNumber == 14) then
open(12, file='Simulated_Data/Simulated_data14.txt')
else if (jExperimentNumber == 15) then
open(12, file='Simulated_Data/Simulated_data15.txt')
else if (jExperimentNumber == 16) then
open(12, file='Simulated_Data/Simulated_data16.txt')
else if (jExperimentNumber == 17) then
open(12, file='Simulated_Data/Simulated_data17.txt')
else if (jExperimentNumber == 18) then
open(12, file='Simulated_Data/Simulated_data18.txt')
else if (jExperimentNumber == 19) then
open(12, file='Simulated_Data/Simulated_data19.txt')
else if (jExperimentNumber == 20) then
open(12, file='Simulated_Data/Simulated_data20.txt')
else if (jExperimentNumber == 21) then
open(12, file='Simulated_Data/Simulated_data21.txt')
else if (jExperimentNumber == 22) then
open(12, file='Simulated_Data/Simulated_data22.txt')
else if (jExperimentNumber == 23) then
open(12, file='Simulated_Data/Simulated_data23.txt')
else if (jExperimentNumber == 24) then
open(12, file='Simulated_Data/Simulated_data24.txt')
else if (jExperimentNumber == 25) then
open(12, file='Simulated_Data/Simulated_data25.txt')
else if (jExperimentNumber == 26) then
open(12, file='Simulated_Data/Simulated_data26.txt')
else if (jExperimentNumber == 27) then
open(12, file='Simulated_Data/Simulated_data27.txt')
else if (jExperimentNumber == 28) then
open(12, file='Simulated_Data/Simulated_data28.txt')
else if (jExperimentNumber == 29) then
open(12, file='Simulated_Data/Simulated_data29.txt')
else if (jExperimentNumber == 30) then
open(12, file='Simulated_Data/Simulated_data30.txt')
else if (jExperimentNumber == 31) then
open(12, file='Simulated_Data/Simulated_data31.txt')
end if 




do jfix=1, min(n_fix, I_state_fix) ; do je = 1, n_e ; do jht=1, min(n_ht, I_state_ht)
    call TYPE1_DET(je,jht,jfix,jtype1)
do jsample=1,int(n_samp(jtype1))   
    
    !read same random draws 
            do jno=1, 443
            read (14,*) ran_hypo(jno) ! 
            end do
                        do jno=1, 213 ! for old
                        read (15,*) ran_hypo_old(jno) ! 
                        end do
                        
    
do jage=1,n_age 

! *************   FIRST PERIOD WHEN AGE=1    **************
if (jage.eq.1) then 
               !jas=1
               jassets=assets_initial(je)
               jx=0 ! 0 years of experience
               jx_int=0.d0    
!Pick jh from the initial distribution of health at age 25 from data
                    rb=ran_hypo(1)
                    bb=1
                    prob=In_h(je,bb,jht)
                 57 if (rb.lt.prob) then
                    jh=bb
                    else 
                    bb=bb+1
                    if (bb.le.n_h) then
		            prob=prob+In_h(je,bb,jht)
		            go to 57
		            else 
		            jh=n_h 
                    end if 
                    end if                 
                    
!Pick jr from the initial distribution of health at age 25 from data
                    rb=ran_hypo(2) 
                    bb=1
                    prob=In_r(je,bb)
                 58 if (rb.lt.prob) then
                    jr=bb
                    else 
                    bb=bb+1
                       if (bb.le.n_r) then
		            prob=prob+In_r(je,bb)
		            go to 58
                    else 
		            jr=n_r 
                    end if 
                    end if 
                    
!Pick jmar from the initial distribution at age 25 from data
                    rb=ran_hypo(3)                
                    prob=In_mar(je,jh)
                    if (rb.le.prob) then
                    jmar=2
                    else 
                    jmar=1
                    end if
                    
!draw spouse work status
                    
                    if (jmar.eq.1) then
                        jos=1 !  does not matter
                    else if (jmar.eq.2) then ! draw spouse work status
                         rb=ran_hypo(4) 
                        if (rb.lt.(prob_os(jage,je,jh,1,jfix))) then
                        jos=1
                        else 
                        jos=2
                        end if 
                    end if 
                    
!Draw employment offer
                  rb=ran_hypo(5) 
		            bb=1
		            prob=Pr_o_25(bb, je)
85		            if (rb.lt.prob) then
		            jo=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_o) then
		            prob=prob+Pr_o_25(bb, je)
		            go to 85
                    else 
                        jo=n_o
		            end if
                    end if
		            
     
                     rb=ran_hypo(6)  
		            bb=1
		            prob=p_zt(je,bb,2)
86		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 86
                    else 
                        jzt=n_zt
		            end if
                    end if             
                    
end if  !for jage=1         

! *************   ALL AGES<N_RET    **************
if (jage.lt.n_ret) then	
    
    !Draw mortality shock            
                     rb=ran_hypo(7+11*(jage-1))
		             if (rb.gt.p_surv(jh,je,jmar,  jage)) then !individual dies	       
                     go to 33 !simulate new individual 
                     end if

call TYPE4_DET(jmar,jos,jtype4) ! if jage is not 1, then jmar and jos have been already drawn towards the end of the code.                                            
call TYPE3_DET(jh,jr,jtype3) ! we need this index only for non-retirees

! determine optimal employment status "jo_act"
   !The individual decides whether to accept the offer or decline
   !interpolate with respect to jas and jx.

 if (jo.eq.1) then 
         jo_act=1  ! THIS IS YOUR ACTUAL EMPLOYMENT STATUS 
         
 else  ! HAVE TO DECIDE BETWEEN WORKING AND NOT WORKING
         EXPECTED_B(jo)=0.d0
         EXPECTED_B(1)=0.0d0

call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,   B(1,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(1))   
 
call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,  B(jo,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(jo))              



         If   (EXPECTED_B(jo).gt.EXPECTED_B(1)) then  
           jo_act=jo   
         else 
           jo_act=1 
         end if 
                
 end if 
 
!Draw health shocks (TYPE2)
                   rb=ran_hypo(8+11*(jage-1))        
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         90 if (rb.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 90
		            else 
		            jtype2=n_type2 
		            end if
                    end if
                    
if (jtype2.eq.1) then
jdu=1 ; jdp=1 ; js=1
else if  (jtype2.eq.2) then
jdu=1 ; jdp=2 ; js=1
else if  (jtype2.eq.3) then
jdu=1 ; jdp=1 ; js=2
else if  (jtype2.eq.4) then
jdu=1 ; jdp=2 ; js=2
else if  (jtype2.eq.5) then
jdu=2 ; jdp=1 ; js=1
else if  (jtype2.eq.6) then
jdu=2 ; jdp=2 ; js=1
else if  (jtype2.eq.7) then
jdu=2 ; jdp=1 ; js=2
else if  (jtype2.eq.8) then
jdu=2 ; jdp=2 ; js=2
end if
		             
!draw the medical expenditure catastrophic shock
                   rb=ran_hypo(9+11*(jage-1)) 
		            bb=1
		            prob= Prob_MEcat(bb)
		         28 if (rb.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 28
		            else 
		            jMEcat= n_MEcat  
		            end if
                    end if  
                    
                    
                    ! draw whether you have the option to not pay and treat. 
                    if ((jo_act.eq.3).OR.(jo_act.eq.5)) then

                        if ((jExperimentNumber.eq.10).OR.(jExperimentNumber.eq.30).OR.(jExperimentNumber.eq.8)) then
                          joption_pay= 1 
                        else 
                           joption_pay=2 ! for those with ESHI, they all have option not to pay and get treated                         
                        end if 
                        
                    else 
                        
                   rb=ran_hypo(10+11*(jage-1)) 
            		prob= Shock_type(1,jo_act,jh,jage) + Shock_type(2,jo_act,jh,jage)   
                  if (rb.le.prob) then
		            joption_pay=1 ! cannot treat and not pay
		            else 
		            joption_pay=2 ! can treat and not pay
                    end if  
                    end if 
                    
                    
                    if ((joption_pay.eq.1).and.(jtype2.ne.1)) then
                    rb=ran_hypo(11+11*(jage-1))
                    if (rb.lt.(Shock_type(1,jo_act,jh,jage)/(Shock_type(1,jo_act,jh,jage)+Shock_type(2,jo_act,jh,jage)))) then ! you cannot treat 
                        juntreated=1 ! cannot treat
                    else 
                        juntreated=2 ! this means you have option to treat
                    end if 
                    else 
                        juntreated=2 ! set to 2 for everyone else 
                    end if 
                    
         if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
 

    if (jo_act.eq.1) then 
     jx_int_p=jx_int
     else if (jo_act.eq.2.OR.jo_act.eq.3) then  
     jx_int_p= min((jx_int+ (max( (hrs_pt - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) )
     else if (jo_act.eq.4.OR.jo_act.eq.5) then 
     jx_int_p=   min((jx_int+ (max((hrs_ft - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) ) ! jx_int goes up by 1 if you work full time and have no sick days. We cap it at 100. 
     end if
     
if (jage.le.n_ret-2) then    
do jzhc_p=1,n_zhc 
      HK_INTERPOL_sim(jzhc_p) =max(0.0d0,  min (jx_int_p *zhc(jzhc_p), dble(n_hc))) 
end do       
end if

Z_INT_OFFER=0.0d0
Z_INT=0.0d0

call InterpolateScalar(HK(0 : n_hc), z(0 : n_hc,jh,je,jfix,jzt), jx_int, Z_INT_OFFER)
    if (jo_act.ge.2) then  
     Z_INT=Z_INT_OFFER
    end if
    
      if (((jExperimentNumber.eq.25).AND.(jage.eq.6)).OR.((jExperimentNumber.eq.26).AND.(jage.eq.16)).OR.((jExperimentNumber.eq.27).AND.(jage.eq.26)).OR.((jExperimentNumber.eq.28).AND.(jage.eq.36))) then
       Z_INT= Z_INT*1.02d0
       Z_INT_OFFER= Z_INT_OFFER*1.02d0
     end if 
   
    
if (jage.eq.(n_ret-1)) then 
!Determine human capital group for the purpose of SS next period for those aged 64 (jhc_ret)
if (jx_int.le.hc_TS(1,je)) then 
    jhc_ret=1
else if (jx_int.le.hc_TS(2,je)) then 
    jhc_ret=2
else 
    jhc_ret=3
end if
else 
   jhc_ret=1  ! just to have a value although this will not matter for anything. 
end if 
 

! calculate incomes if pay or not pay
 actual_earn=0.0d0
 income_bt=0.0d0
 tax_SS=0.0d0
 tax_med=0.0d0
 
  actual_earn = max(0.0d0 , ( wage_penalty(jo_act,je)*Z_INT * (hrs_offer(jo_act) -phi(je,jh,jtype2) )) ) ! earnings net of sick days
  income_bt= ((interest-1.0d0))*jassets +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
 
    if ((jExperimentNumber.eq.14).OR.(jExperimentNumber.eq.15).OR.(jExperimentNumber.eq.23)) THEN
      
        if (jo.ne.1) then ! husband works - deduct from his income
            tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo_act,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo_act,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife.
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)))
end if 
            
            else 
if ((jo_act.eq.3).OR.(jo_act.eq.5)) then ! husband holds insurance -- so deduct from his income only
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo_act,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo_act,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. note this can be zero
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)))
end if   

end if 
  
 ! for all of these, pay OOP =1, not pay OOP=2 
 income_pay(:)=0.0d0
 Total_tax_pay(:)=0.0d0
 base_2_pay(:)=0.0d0
 base_1_pay(:)=0.0d0
 
        ! if pay OOP 
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo_act) + oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo_act,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo_act) +  oop_spouse(jage,je,jtype4,jo_act) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)
      ! if not pay OOP
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(2) =  prem_ephi(jo_act,jtype4) +  oop_spouse(jage,je,jtype4,jo_act) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)
  
  
  
!Determine income group
          if (income_bt.le.INC_TS(1)) then 
        jinc=1
        else if (income_bt.le.INC_TS(2)) then 
        jinc=2
        else if (income_bt.le.INC_TS(3)) then 
        jinc=3
        else if (income_bt.le.INC_TS(4)) then 
        jinc=4
        else 
        jinc=5
        end if
    

        
!DETERMINE IF THE PERSON WILL PAY/TREAT and OPTIMAL CONSUMPTION/SAV

  !initiate        
CON_sim(:)=0.0d0
  ASP_sim(:)=0.d0
  B_sim(:)=0.d0
  TR_sim(:)=0.0d0
  
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay=0.0d0
  V_c_all_trpay=0.0d0
  V_death=0.0d0
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
        UTI_print=0.0d0
        V_death_print=0.0d0
  
!If I don't have a shock, I pay medical bill. 
if (jtype2.eq.1) then
     jtrpay=1
    I_treat=2
    I_pay=1
    
    !If on cons floor
    if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then 
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    
      ! still need to calculate V because we need it all the time for welfare analysis
    ! get utility
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME   
    else
              jtr_act=0.d0 
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if
               
!VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets   
                    ! the value function tomorrow depends on whether I'm working today
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,1,je,jh,jo,jtype4,jtype2,(base_1_pay(1)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 

          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit  
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,1,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               end if 
               
                  EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 

     if (B_act.lt.maxB) exit
              maxB=B_act 
              jas_p_act=ASP
              jcons_act=CON
    
       
51 end do !jsav
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  
      
           B_sim(1) = maxB
            
    end if 
    
    
    
 !If I have a shock   
else 
if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.(juntreated.EQ.1).AND.(I_treat_Mcaid.EQ.0)) then ! if I'm on the cons floor even if I don't pay AND cannot treat
    

        jtrpay=3 ! you don't treat
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(2)
  

    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 3
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,1,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(2,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
    
    
else if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.((juntreated.EQ.2).OR.(I_treat_Mcaid.EQ.1))) then ! if I'm on the cons floor even if I don't pay AND can treat - then you always "treat and pay"
    jtrpay=1
    I_treat=2
    I_pay=1
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    

    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
else
    
if (joption_pay.eq.2) then   ! if I have all 3 options 
do jtrpay=1, 3
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else               
          EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      end if 
      
     end do ; end do ; end do   ; end do; end do  
  
         end if 
         
         
         
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    
            
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets  
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
      if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat              
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else
                              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

      end if 
      
      
      end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 
 
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.   
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     
     else
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit  
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                                        if (jo_act.eq.1) then

                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                                        else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                        end if 
                                        
                    
                 
                    
               if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat     
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
               else 
                                       EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

               end if 
               
               
               end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
  
end do ! jtrpay

! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).OR.(B_sim(2).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 


if ((B_sim(1).gt.B_sim(2)).AND.(B_sim(1).ge.B_sim(3))) then
    jtrpay = 1 
else if ((B_sim(2).ge.B_sim(1)).AND.(B_sim(2).ge.B_sim(3))) then
    jtrpay = 2 
else if ((B_sim(2).eq.B_sim(1)).AND.(B_sim(2).eq.B_sim(3))) then
    jtrpay = 2 
    else
    jtrpay = 3 
    end if 
    
 if (jExperimentNumber.eq.13) THEN ! here we force people to pay and treat
      jtrpay = 1 
 end if 
 
    

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)

 else if (((joption_pay.eq.1).AND.(juntreated.EQ.2)).OR.(((joption_pay.eq.1).AND.(juntreated.EQ.1)).AND.(I_treat_Mcaid.EQ.1).AND.(base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) )) then ! if you don't have the option to treat and not pay. BUT HAVE OPTION TO TREAT

    do jtrpay=1, 3, 2 ! i only consider 1 and 3
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then 
 jcons_act=cons_floor(jage, jmar,je,jh)  
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    
               
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets    
                     if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                     else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                             end if 
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    end if 
                    
                        
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav

   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
    end do ! jtrpay
    
    
    ! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 
if (B_sim(1).ge.B_sim(3)) then
    jtrpay = 1 
    else
    jtrpay = 3 
end if 

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)


 else if ((joption_pay.eq.1).AND.(juntreated.EQ.1)) then ! you cannot treat
    
    
    jtrpay=3 ! you don't treat
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    

        I_treat=1 ! not treat
        I_pay=2 ! not pay


!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then 
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    
               
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
                             call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
                    
                        
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
  
  
jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay) 
    
    

    
    
end if ! for having the options in trpay


    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else if (jtrpay.eq.3) then
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

    

    
end if ! for being on consumption floor even if you don't pay
end if ! for having a health shock

!If you don't get treated, then there is a probability the shock and R don't get reported

!initialize to not having shock
jdp_new=jdp
jdu_new=jdu
js_new=js
jr_rep=jr

!dp (2,4,6,8)
!if ((jtype2.eq.2).OR.(jtype2.eq.4).OR.(jtype2.eq.6).OR.(jtype2.eq.8) ) then
if (jdp.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,1,I_treat))   then
      jdp_new=1 
    end if 
end if 


!du (5,6,7,8)
!if ((jtype2.eq.5).OR.(jtype2.eq.6).OR.(jtype2.eq.7).OR.(jtype2.eq.8) ) then
 if (jdu.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,2,I_treat))   then
      jdu_new=1 
    end if 
end if

!s (3,4,7,8)
 if (js.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,3,I_treat))   then
      js_new=1 
    end if 
 end if
 
 ! R
 if (jr.eq.2) then
         rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_R(jo_act))   then
      jr_rep=1 
    end if 
 end if 
 

 if (jtr_act.gt.0.0d0) then 
 jdp_new=jdp
 jdu_new=jdu
 js_new=js
 jr_rep=jr
 end if 
 

call TYPE2_DET(jdu_new,jdp_new,js_new,jtype2_new)


if (jtype2.eq.1) then
         call func_uti(jage,1,je,jh,jo_act,jtype4,jtype2,jcons_act,UTI_print) 
else 
         call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,jcons_act,UTI_print) 
end if 

        V_death_print =  Bequest_val(jas_p_act)  + cost_death


!Write to file
 write(12,1008) jsample, jage, je, jfix, jht, jmar, jh,  jr, jtype2, jtype4, jo ,   jo_act,  jzt, jx_int, jx_int_p, oop(jtype2, jh, jage,jMEcat,jo_act)*5200.0d0 , Charges(jtype2, jh, jage,jMEcat)*5200.0d0, jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , wage_penalty(jo,je)*Z_INT_OFFER   , 5200.0d0*actual_earn, max(0.0d0, (hrs_offer(jo_act) -phi(je,jh,jtype2) ))*100.0d0, exp( log(wage_penalty(jo_act,je)*Z_INT_OFFER)+random_normal()*stdze(je) ) ,  inc_spouse(jage,je,jh,jtype4,jfix)*5200.0d0, jtrpay  , oop_spouse(jage,je,jtype4,jo_act)*5200.0d0 ,base_1_pay(1)*5200, base_1_pay(2)*5200, base_2_pay(1)*5200, base_2_pay(2)*5200, 0, jtype2_new, jr_rep, joption_pay, juntreated, B_sim(jtrpay), UTI_print, V_death_print
 
    
!draw marital status for next period 
            rf=ran_hypo(12+11*(jage-1)) 
            if (rf.lt.(Tp_mar(1,jmar,je,jage,jh,jinc,jo_act))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p   
!Health Risk R
          rb=ran_hypo(13+11*(jage-1)) 
          if (rb.lt.p_risk(jr, 1,jage,jh)) then
            jr_p=1
            else 
            jr_p=2
            end if
          jr=jr_p    
!Health
         rb=ran_hypo(14+11*(jage-1)) 
            bb=1
            prob=Tp_h(bb, jh,je,jht,jage, jtype2,I_treat,jmc) 
         92 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
                if (bb.le.n_h) then
                prob=prob+Tp_h(bb, jh,je,jht,jage, jtype2,I_treat,jmc)
                go to 92
                else 
              jh_p= n_h  
                end if
            end if  
            jh=jh_p        
            
!Next period states needed only for those of working age, but not 64: 
if (jage.NE.(n_ret-1)) then    
    
    
                
        jx_int=max(0.0d0, min(jx_int_p , dble(n_hc)) )
            
!Draw employment offer jo

     ro=ran_hypo(15+11*(jage-1)) ! draw a random # between 0-1
        bb=1
        prob=Pr_o_i(bb,jo_act, je,jage)   
     34 if (ro.lt.prob) then
        jo=bb
        else 
        bb=bb+1
        if (bb.le.n_o) then
        prob=prob+Pr_o_i(bb,jo_act, je,jage)
            
        go to 34
        else 
            jo=n_o
        end if
        end if   
        
        
        
!Draw zt next period
         rb=ran_hypo(16+11*(jage-1)) ! draw a random # between 0-1
        if (jo_act.eq.1) then                    
		            bb=1
		            prob=p_zt(je,bb,1)
89		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,1)
		            go to 89
                    else 
                        jzt=n_zt
		            end if
                    end if  
  
        else 
 
                            
		            bb=1
		            prob=p_zt(je,bb,2)
99		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 99
                    else 
                        jzt=n_zt
		            end if
                    end if  
        end if 
        
!Spouse's work status next period (jos)- this depends on the individual's state next period, so it's ok to enter the new jh here. Also, we've already determined the new jmar next period

                        ! draw spouse work status  1=not work, 2=work ; if not married, doesn't matter, jos will be drawn but not used.
                        rb=ran_hypo(17+11*(jage-1))
                        if (rb.lt.(prob_os(jage+1,je,jh,1,jfix))) then ! this is the jh next period   
                        jos=1
                        else 
                        jos=2
                        end if 
               
end if ! for age not being 64

end if ! for jage<n_ret       


		
!****************************************************************************************************************************************


if (jage.ge.n_ret) then !individual is retired	
           
!Draw health shocks jtype2
                  rf=ran_hypo_old((jage-41)*6+1)
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         37 if (rf.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 37
                    else 
                        jtype2=n_type2
		            end if
                    end if		         
!Draw survival		            
                    rx=ran_hypo_old((jage-41)*6+2)
		             if (rx.gt.p_surv(jh,je,jmar , jage)) then !individual dies	
                     go to 33 !simulate new individual 
                     end if	             
!Draw the medical expenditure catastrophic shock
                    rf=ran_hypo_old((jage-41)*6+3)
		            bb=1
		            prob= Prob_MEcat(bb)
		         18 if (rf.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 18
                    else 
                        jMEcat=n_MEcat
		            end if
                    end if            
                     
		             


!income
income = max(0.0d0, (interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -  prem_mcare(jmar)  - max(0.d0,  oop(jtype2, jh, jage,jMEcat,1) +  oop_spouse(jage,je,jmar,1) - 0.075d0*((interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )) !taxable income. the last term is OOP in excess of a threshold of your income, which is deductable.
base_2 = prem_mcare(jmar) + oop(jtype2, jh, jage,jMEcat,1) + oop_spouse(jage,je,jmar,1)  + max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) ! total resources subtracted: premiums, OOP, and tax, at the HH level

                
        
!Initialize
 RET_act=0.0d0
 jtr_act=0.0d0
 jcons_act=0.0d0
 jas_p_act=0.0d0
 maxB=0.0d0
 UTI_print=0.0d0
 V_death_print=0.0d0
  
if (((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) .le. (cons_floor_at(jage, jmar,je,jh)) ) then    ! get transfer                                                                                                                                  
    jtr_act= cons_floor_at(jage, jmar,je,jh) - (interest*jassets+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2) 
    jcons_act=cons_floor(jage, jmar,je,jh)
    jas_p_act=0.0d0

else
   jtr_act=0.d0
!Find max savings point
   MAX_SAV=  interest*jassets  + pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2  -cons_floor_at(jage, jmar,je,jh) -jassets
   
!Find lower bound
            jsav_upper_ret=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper_ret+1) )
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper_ret.lt.1) then
                 PRINT *,'Error'
               end if 
               
               
               
!Find value of consuming all and saving nothing
   CON=((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)      
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
                        V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         
                if (jage.eq.n_age) then
                 V_c_all_trpay=UTI+ beta_hat(je)* V_death  
                else 
                   EXPECTED_RET_PRIME=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do 
                 
                 V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
                end if 
         
!Find value of consuming the consumption floor and saving as much as possible 
      call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
      ASP= MAX_SAV+jassets
      V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death    
                   
                if (jage.eq.n_age) then
                 V_cf_trpay=UTI+ beta_hat(je)* V_death  
                else 
                  EXPECTED_RET_PRIME=0.0d0   
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        V_cf_trpay=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                 
               

!loop for optimal savings
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
   do jsav= jsav_upper_ret , 1, -1
       
!Find CON, ASP and UTI

    ASP=  jassets+  sav_sim(jsav)
                if (ASP.lt.0.0d0)  exit
                V_death=cost_death + Bequest_val(ASP)
                
        CON= ( ((interest-1.0d0))*jassets  + pen(jhc_ret,  je,jfix)+ inc_spouse(jage,je,1,jmar,jfix) -base_2 -sav_sim(jsav) )/(1.0d0+taxc)

        if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
         PRINT *,'Error'
          pause
        end if      
 
     call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
     
!Find continuation value
            
                EXPECTED_RET_PRIME=0.0d0 
                
                if (jage.eq.n_age) then
                 RET_act=UTI+ beta_hat(je)* V_death  
                else 
                    
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        RET_act=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                
!Find max value
                if (RET_act.lt.maxB) exit
                maxB=RET_act 
                jas_p_act=ASP
                jcons_act=CON         
    		    	    
56 end do !for jsav for retired ppl
   
     if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= ((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)
     end if 
     
end if ! for being on tr
      

call func_uti_ret(je,jh,jmar,jtype2,jcons_act,UTI_print)
V_death_print=cost_death + Bequest_val(jas_p_act)
   
!Write decisions to file
write(12,1008)  jsample, jage, je, jfix, jht, jmar, jh,  jr, jtype2, 0, 0 ,  0,  0, 0.0d0 , 0.0d0, oop(jtype2, jh, jage,jMEcat,1)*5200.0d0 , 0.0d0 , jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , 0.0d0   , pen(jhc_ret,  je,jfix)*5200.0d0 ,  0.0d0 ,  0.0d0,  inc_spouse(jage,je,1,jmar,jfix)*5200.0d0, 0  , oop_spouse(jage,je,jmar,1)*5200.0d0, ((interest*jassets+  pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 ))*5200.0d0 , 0.0, base_2*5200.0d0, 0.0, jhc_ret, 0, 0, 0, 0, maxB, UTI_print, V_death_print

if (jage.eq.n_age) then
go to 33 !simulate new individual      
end if                            
      
! draw marital status for next period
            rf=ran_hypo_old((jage-41)*6+4)
            if (rf.lt.(Tp_mar(1,jmar,je,jage,1,1,1))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p            
!Risk factors
          rb=ran_hypo_old((jage-41)*6+5)
            bb=1
            prob=p_risk(jr, bb,jage,jh)
         39 if (rb.lt.prob) then
            jr_p=bb
            else 
            bb=bb+1
            if (bb.LE.n_r) then
            prob=prob+p_risk(jr, bb,jage,jh)
            go to 39
            else 
            jr_p=n_r  
            end if
            end if
           jr=jr_p        
!Health
            rb=ran_hypo_old((jage-41)*6+6)
            bb=1
            prob=Tp_h(bb, jh,je,jht,jage, jtype2,2,1) ! YOU ALWAYS GET TREATED IF RETIRED
         32 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
            if (bb.le.n_h) then
            prob=prob+Tp_h(bb, jh,je,jht,jage, jtype2,2,1)
            go to 32
            else 
            jh_p= n_h 
            end if
            end if 
    	    jh=jh_p          
end if !for being retired
           
!Assets next period  
  if (jage.lt.n_age)then  
            jassets=jas_p_act          
  end if                 		            
  
end do !age
33 end do !sample
   end do; end do; end do ! jht, jfix, je

close (12)
close (14)
close (15)
End Subroutine